\documentclass[11pt]{beamer}
\usepackage[UTF8,scheme=plain]{ctex}
\usepackage{listings}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{graphicx}
\usetheme{Boadilla}

\usepackage{framed} % 可以用 \begin{shaded}，即背景色块
\definecolor{shadecolor}{rgb}{0.9,0.9,0.9}

\lstset{
	basicstyle=\footnotesize,
    flexiblecolumns,
    showstringspaces = false,
	numbers=left,
    numberstyle=\tiny,
    numbersep=4pt,
	breaklines,%自动换行
	frame=shadowbox,%加边框
	backgroundcolor = \color{red!30!green!30!blue!30},
	%背景色：蓝灰色
	language=C++, % 语言为 c++
	commentstyle=\color{red!80!green!80!blue!80},% 注释用蓝灰色字体标出
	keywordstyle=\color{blue!70}, % 高亮部分用蓝色字体
	tabsize = 4,
    escapeinside=@@
}

\begin{document}
	\author{ 路毅 \hspace{0.3cm} 曲阜师范大学 }
	\date{\number\year 年 \number\month 月 \number\day 日}
	\title{计算物理第5部分：蒙特卡洛方法}
	%\subtitle{}
	%\logo{}
	%\institute{}
	%\date{}
	%\subject{}
	%\setbeamercovered{transparent}
	%\setbeamertemplate{navigation symbols}{}
	\begin{frame}[plain]
	\maketitle
\end{frame}

\kaishu

\begin{frame}{蒙特卡洛方法}
欧洲的摩纳哥公国，有一个美丽的城市叫做蒙特卡洛。
\begin{figure}
\centering
\includegraphics[width=0.7\textwidth]{MonteCarlo.jpeg}
\end{figure}
\end{frame}

\begin{frame}{蒙特卡洛方法}
蒙特卡洛是著名的赌城，每天都进行大量的赌博。
\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{赌圣.jpeg}
\end{figure}
\end{frame}

\begin{frame}{蒙特卡洛方法}
随机性（运气），是赌博的激情的来源。。。
\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{dice.png}
\end{figure}
\end{frame}

\begin{frame}{蒙特卡洛方法}
\begin{figure}
\centering
\includegraphics[width=0.35\textwidth]{VonNeumann.jpeg}
\end{figure}
冯诺依曼：把使用随机数的一类方法叫做{\bf “蒙特卡洛方法”}。
\end{frame}

\begin{frame}
	\begin{itemize}
		\item {\color{blue} 第十二章：蒙特卡洛方法 }
            \begin{itemize}
            \item {\color{blue} 简介}
            \item 蒙特卡洛模拟示例
            \item 蒙特卡洛积分
            \item Markov链
            \end{itemize}
		\item 第十三章：变分蒙特卡洛
	\end{itemize}
\end{frame}

\begin{frame}[fragile]{随机数的产生，乘加同余法}
传统的计算机不产生随机性。\\
\vspace{0.3cm}
伪随机数：用计算机产生一个无限长的数列，让数列上的数彼此间没有关联，整体均匀分布在某个区间。\\
\vspace{0.3cm}
乘加同余法，
\begin{eqnarray}
I_{j+1} = a I_j + c (\rm mod~~ m),
\end{eqnarray}
其中 $a,c,m$ 均为正整数，产生均匀分布的随机整数，当$m$为很大的整数时具有较大的周期。\\
\vspace{0.3cm}
{\color{blue} 为了方便，从现在开始，我们说的“随机数”都是指“伪随机数”。}
\end{frame}

\begin{frame}[fragile]{随机数的产生，乘加同余法}
在c++中，cstdlib库中的rand()函数产生[0,RAND\_MAX]区间内的随机数，RAND\_MAX的值是一个很大的整数，与所用的机器有关。\\
下面的代码输出RAND\_MAX的值，并产生10个随机数。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cstdlib>

int main(){

        int i;
        cout<<"RAND_MAX="<<RAND_MAX<<endl;
        for(i=0;i<10;i++)
                cout<<rand()<<endl;
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{随机数的产生，乘加同余法}
srand(time(NULL))：将随机数种子设定为1970年1月1日0点到运行程序的那一刻的总秒数，
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cstdlib>

int main(){

        srand(time(NULL)); // 随机数种子：1970年1月1日0点到此刻的总秒数
        int i;
        cout<<"RAND_MAX="<<RAND_MAX<<endl;
        for(i=0;i<10;i++)
                cout<<rand()<<endl;
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{离散的随机数}
rand()对15求余数，得到在$0,1,\cdots,14$中的离散随机数，
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<string>
#include<time.h>
#include<vector>

const string name[15]={"邓琦","李雪","刘浩","吕彦萌","王庆轩","王群","王小鹏","王妍","杨浩","杨红艳","张晓璐","赵翔宇","甄洪悦","周琪琪","张莹"};
int main(){

        int i,j,n=-1;
        srand(time(NULL));
        cout<<"你要点多少次名？ n="; cin>>n;
        for(j=0;j<n;j++){
                i = rand() % 15;
                cout<<j+1<<"\t"<<name[i]<<endl;
        }
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{“连续”分布的随机数}
我们也可以构造$[0,1]$区间内随机分布的“连续”随机数，
\begin{lstlisting}
(double)rand()/RAND_MAX
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{产生任意概率分布的随机数}
构造随机数：$x\in[a,b]$，其概率密度函数为$g(x)$，
\begin{equation}
\int^b_a g(x) dx = 1,
\end{equation}
\begin{shaded*}
\begin{itemize}
\item [1] 找到一个值$Max$，使得$\forall x \in [a,b]$，$g(x) \leq Max$。
\item [2] 构造随机数$\xi$，在$[0, Max]$均匀分布
            \lstset{numbers=none}
            \begin{lstlisting}
            (double)rand() / RAND_MAX * Max
            \end{lstlisting}
\item [3] 构造随机数$x$，在$[a,b]$均匀分布
            \lstset{numbers=none}
            \begin{lstlisting}
            (double)rand() / RAND_MAX * (b-a) + a
            \end{lstlisting}
\item [4] 每取一个$x$，也取一个$\xi$；如果$\xi < g(x)$就使用这个$x$；如果$\xi > g(x)$就不使用这个$x$。
\end{itemize}
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{产生任意概率分布的随机数}
\begin{lstlisting}
/*
 * random(...)  产生在 [a,b] 区间内，以 func(x) 为概率密度的随机数，在使用这个函数之前，可以用srand()
函数对随机数种子进行初始化
 * double a
 * double b                     随机数所在区间 [a,b]
 * double (*func)(double x)     随机数的相对概率密度
 * double fmax                  相对概率密度的上界
 */
double gnr_rand(double a, double b, double (*func)(double x), double fmax){        
        double x = (double)rand() / RAND_MAX * (b-a) + a; // 在 [a,b] 之间均匀分布
        double f = (double)rand() / RAND_MAX * fmax; // 在 [0, fmax] 之间均匀分布
        while( f > func(x) ){ // x 有 func(x) 的相对概率被接受
                x = (double)rand() / RAND_MAX * (b-a) + a;
                f = (double)rand() / RAND_MAX * fmax;
        }
        return x;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{检验随机数的“频率”}
另外我们还可以写一个函数，对随机数在各个区间出现的频率进行统计
\begin{lstlisting}
/*
 * count_random(...) 如果rand出现在 [a+i*h, a+(i+1)*h)中，那么count[i]加1
 */
void count_random(double a, double h, int * count, double rand){
        
        count[ (int)((rand-a)/h) ] ++;
}  
\end{lstlisting}
我们可以产生$[-3,3]$区间内服从标准高斯分布$N(\mu=0, \sigma^2=1)$的随机数
\begin{equation}
f(x) = \frac{1}{\sqrt{2\pi}} exp(-\frac{x^2}{2}),
\end{equation}
\end{frame}

\begin{frame}[fragile]{检验随机数的“频率”}
并且统计产生的随机数，与概率密度函数相比照，检验我们的代码
\begin{lstlisting}
头文件。。。

double f(double x){//标准高斯分布，或标准正态分布
        return exp(-x*x/2)/sqrt(2*M_PI);
}

int main(){

        int i,N=1000000;
        double a=-3, b=3, fmax = 1/sqrt(2*M_PI);
        double h=0.1, x;
        int n = (b-a)/h;
        int *count = new int [n];
        for(i=0;i<n;i++) count[i] = 0;
        for(i=0;i<N;i++){
                x = gnr_rand(a, b, f, fmax);
                count_random(a, h, count, x);
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{检验随机数的“频率”}
\lstset{firstnumber=last}
\begin{lstlisting}
        //输出统计“频率”
        ofstream fp("Gauss_rand.txt");
        for(i=0;i<n;i++){
                fp<< a+ (i+0.5)*h <<"\t"<<(double)count[i]/N/h<<endl;
                //第i个区间内的数的频率除以区间长度，得到“频率密度”
        }
        fp.close();
        delete [] count;
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{检验随机数的“频率”}
得到“频率密度”数据，以后，可以与概率密度函数相比较，gnuplot代码为
\begin{lstlisting}
plot exp(-x*x/2)/sqrt(2*3.1415926), "Gauss_rand.txt" w p
\end{lstlisting}
得到下图，可以看到，我们产生的随机数很好地符合高斯分布。
\begin{figure}
\centering
\includegraphics[width=0.6\textwidth]{Gauss_rand}
%\caption{标准高斯分布，随机数“频率密度”与概率密度的比照}
\end{figure}
如果需要其他分布的随机数，只需修改程序中的f(x)和fmax即可。
\end{frame}

\begin{frame}[fragile]{大数定理与中心极限定理}
随机数$X$服从任意一个分布，其概率密度函数为$f(x)$，期望值为
\begin{equation}
E(X) = \int^\infty_{-\infty} x f(x) dx,
\end{equation}
方差为
\begin{equation}
D(X) = \int^\infty_{-\infty} (x - E(X))^2 f(x) dx = \sigma^2.
\end{equation}
那么，如果每次抽样取$n$个独立服从$f(x)$的随机数，然后取其平均数，则有
\begin{eqnarray}
\lim_{n \rightarrow \infty} P \left\{ | \frac{1}{n} \sum\limits^n_{i=1} X_i - \mu | \leq \epsilon \right\}
\geq \lim_{n \rightarrow \infty} \left[1- \frac{ \sigma^2 }{ n \epsilon^2} \right]
 =1.
\end{eqnarray}
这就是{\bf 大数定理},它的意思很直白：随着抽样次数的增大，抽样平均一定会收敛于期望值。
方差越小，收敛得越快。
\end{frame}

\begin{frame}[fragile]{大数定理与中心极限定理}
中心极限定理则不那么直白，但是可能更神奇：
\begin{eqnarray}
\lim_{n \rightarrow \infty} \frac{ \frac{1}{n}\sum\limits^n_{i=1} X_i - \mu }{ \sigma / \sqrt{n} } \sim N(0,1),
\end{eqnarray}
它的意思是说：如果抽样次数足够大，抽样平均服从正态分布$N(\mu, \sigma^2/n)$。
\begin{itemize}
\item[1] 不论$X$服从什么分布，不论$f(x)$是什么样的概率密度分布函数，抽样次数足够多时$X$的抽样平均都服从正态分布。
\item[2] 所以正态分布在大自然中具有特殊的意义，它是分布中的分布。
\end{itemize}
这两个定理在后文中我们会用到，现在我们不妨用刚学会的伪随机数，稍微检验这两个定理。
\end{frame}

\begin{frame}[fragile]{检验中心极限定理}
我们可以使用$[-3,3]$之间的均匀分布，做30次抽样，取平均，然后统计它的“频率”密度，再与$N(\mu, \sigma^2/30) = N(0,0.1)$做比较，从比较图上便可以一目了然地看到，符合得好不好。
下面的代码实现这个比较
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cstdlib>
#include<cmath>
#include<fstream>

/*
 * random(...)  产生在 [a,b] 区间内，以 func(x) 为概率密度的随机数，在使用这个函数之前，可以用srand()函数对随机数种子进行初始化
 * double a
 * double b                     随机数所在区间 [a,b]
 * double (*func)(double x)     随机数的相对概率密度
 * double fmax                  相对概率密度的上界
 */
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{检验中心极限定理}
\lstset{firstnumber=last}
\begin{lstlisting}
double gnr_rand(double a, double b, double (*func)(double x), double fmax){

        double x = (double)rand() / RAND_MAX * (b-a) + a; // 在 [a,b] 之间均匀分布
        double f = (double)rand() / RAND_MAX * fmax; // 在 [0, fmax] 之间均匀分布
        while( f > func(x) ){ // x 有 func(x) 的相对概率被接受
                x = (double)rand() / RAND_MAX * (b-a) + a;
                f = (double)rand() / RAND_MAX * fmax;
        }
        return x;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{检验中心极限定理}
\lstset{firstnumber=last}
\begin{lstlisting}
double sample(int n, double a, double b, double (*func)(double x), double fmax){

        double y=0;
        for(int i=0;i<n;i++){
                y += gnr_rand(a, b, func, fmax);
        }
        return y/n;
}

/*
 * count_random(...) 统计随机数rand出现在 [a,a+h), [a+h, a+2h), ... 中的频率，记录在count中
 */
void count_random(double a, double h, int * count, double rand){

        count[ (int)((rand-a)/h) ] ++;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{检验中心极限定理}
\lstset{firstnumber=last}
\begin{lstlisting}
double f(double x){//[-3:3]之间均匀分布
        if(x >= -3 && x <= 3) return 1.0/6;
        else return 0;
}

double NormDistr(double mu, double sigma, double x){
        return 1 / sqrt(2*M_PI) / sigma * exp( - (x-mu)*(x-mu)/2/sigma/sigma );
}

int main(){

        int i,N=1000000;
        double a=-3, b=3, fmax = 1.0/6;
        double h=0.01, x;
        int n = (b-a)/h;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{检验中心极限定理}
\lstset{firstnumber=last}
\begin{lstlisting}
        int *count = new int [n];
        for(i=0;i<n;i++) count[i] = 0;
        for(i=0;i<N;i++){
                x = sample(30, a, b, f, fmax);
                count_random(a, h, count, x);
        }
        ofstream fp("test_LLNandCLT.txt");
        for(i=0;i<n;i++){
                fp<< a+ (i+0.5)*h <<"\t"<<(double)count[i]/N/h<<"\t"<< NormDistr(0, sqrt(3.0/30), a+ (i+0.5)*h )<<endl;
        }
        fp.close();
        delete [] count;

        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{检验中心极限定理}
得到数据文件以后，作图
\begin{lstlisting}
plot "test_LLNandCLT.txt" u 1:2, "" u 1:3 w l
\end{lstlisting}
得到下图
\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{test_LLNandCLT}
\caption{检验大数定理与中心极限定理。用$[-3,3]$之间均匀分布的随机数，进行30次抽样并做平均，根据中心极限定理，这个抽样平均服从$N(0,0.1)$分布。绿色折线为$N(0,0.1)$的概率密度函数，蓝色点为抽样平均值的“频率”密度。}
\label{fig:LLNandCLT}
\end{figure}
\end{frame}

\begin{frame}
	\begin{itemize}
		\item {\color{blue} 第十二章：蒙特卡洛方法 }
            \begin{itemize}
            \item {\color{blue} 简介}
            \item {\color{blue} 蒙特卡洛模拟示例}
            \item 蒙特卡洛积分
            \item Markov链
            \end{itemize}
		\item 第十三章：变分蒙特卡洛
	\end{itemize}
\end{frame}

\begin{frame}[fragile]{随机模拟:气体扩散}
\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{mc_gas}
\caption{蒙特卡洛模拟：气体扩散的简化模型：任一时刻任一气体分子都有一半可能穿过小孔到另一侧。}
\label{fig:Gauss_rand}
\end{figure}
有一个封闭腔，如上图，初始时，有1千个气体分子在隔板左侧，此后，任意时刻随机选取一个气体分子穿过小孔到达另一边，那么一段时间之后，我们可以预期，气体分子在两边的分布会是均匀的。
\end{frame}

\begin{frame}[fragile]{随机模拟:气体扩散}
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>

int main(){

        int i, j, n=1000, Nt = 10000;
        int * walker = new int [n];
        for(i=0;i<n;i++){
                walker[i] = -1;// -1/1 for left/right
        }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{随机模拟:气体扩散}
\lstset{firstnumber=last}
\begin{lstlisting}
        int count[2]={1000,0};
        ofstream fp("mc_gas.txt");
        for(j=1;j<=Nt;j++){
                i = rand() % n ;
                walker[i] *= -1;
                if( walker[i] == 1 ){
                        count[1] ++; count[0] --;
                }
                else{
                        count[0] ++; count[1] --;
                }
                fp<<j<<"\t"<<count[0]<<"\t"<<count[1]<<endl;
        }
        fp.close();
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{随机模拟:气体扩散}
得到数据以后，用gnuplot如下作图，
\begin{lstlisting}
plot "mc_gas.txt" u 1:2 w p t "number of molecules on the left", 500 w l notitle
\end{lstlisting}
得到下图，过一段时间以后，两边气体分子数平衡，此后虽然有些涨落，但基本保持平衡。
\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{mc_gas2}
\caption{气体扩散简单模拟}
\label{fig:mc_gas2}
\end{figure}
\end{frame}

\begin{frame}[fragile]{链式反应：原子弹临界尺寸}
U235的链式反应为
\begin{equation}
n + U235 = n + n + else,
\end{equation}
其中$else$表示其他反应产物。
反应后的2个中子如果再与2个U235反应，就会产生4个中子，如此下去，很快就会触发大量反应。
每次反应都会释放MeV量级的反应（生化反应是eV量级），所以大量这样的核反应会产生巨大的能量，迅速产生很高的温度，将材料汽化电离，释放冲击波和大量辐射，这就是1945年在广岛长崎发生的原子弹爆炸。\\
\vspace{0.3cm}
百度百科：“原子弹爆炸的强烈光波，使成千上万人双目失明；6千多度的高温，把一切都化为灰烬；放射雨使一些人在以后20年中缓慢地走向死亡；冲击波形成的狂风，又把所有的建筑物摧毁殆尽。”
\end{frame}

\begin{frame}[fragile]{链式反应：原子弹临界尺寸}
\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{GuangDaoAtomicBomb}
%\caption{1945年广岛原子弹爆炸图片}
\end{figure}
\end{frame}

\begin{frame}[fragile]{链式反应：原子弹临界尺寸}
但是，如果一次反应释放的2个中子与U235再次反应的概率不大，那么就不会触发链式反应，就不会爆炸。
\\
\vspace{0.3cm}
所以如果U235材料的体积特别小，一个U235衰变以后释放的2个中子都飞出了材料，就不会爆炸。
\\
\vspace{0.3cm}
如果U235材料的体积足够大，一次反应产生的2个中子基本都会遇到另2个U235，那么就会爆炸。
\\
\vspace{0.3cm}
所以，原子弹有一个临界尺度，超过尺度就会爆炸，小于尺度就是安全的。
\end{frame}

\begin{frame}[fragile]{链式反应：原子弹临界尺寸}
我们可以用随机数模拟这个过程。
我们将长度单位选为中子在材料中的平均自由程，所以在这个单位下，平均自由程为1。
\\
\vspace{0.3cm}
为了简化处理，我们假设中子在材料中的自由程（即两次碰撞之间飞行的距离）为$[0,2]$内均匀分布的随机数，那么自然平均自由程是1。
\\
\vspace{0.3cm}
假设有球形原子弹，半径为R。
\\
\vspace{0.3cm}
那么，我们在原子弹体积内均匀取1000个随机点，认为这些点发生了自发裂变，出现了1000个自由中子。
\end{frame}

\begin{frame}[fragile]{链式反应：原子弹临界尺寸}
自由中子的飞行方向随机，朝各个立体角都是相同的概率。
所以$(\theta, \phi)$中$\phi$为$[0,2\pi]$内均匀分布的随机数，而$\theta$要使得
\begin{equation}
d \Omega = \sin \theta d \theta d\phi = d \cos \theta d \phi,
\end{equation}
均匀分布，所以$\cos \theta$是$[-1,1]$内均匀分布的随机数。\\
自由飞行之后的坐标为
\begin{eqnarray}
x^\prime &=& x + r \sin \theta \cos \phi, \\
y^\prime &=& x + r \sin \theta \sin \phi, \\
z^\prime &=& x + r \cos \theta, \\
\end{eqnarray}
终点若在弹体内部，则有
\begin{equation}
(x^\prime)^2 + (y^\prime)^2 + (z^\prime)^2 \leq R^2,
\end{equation}
\end{frame}

\begin{frame}[fragile]{链式反应：原子弹临界尺寸}
我们可以分别计算这1000个中子的自由程；如果落在原子弹内部，就触发一次核反应，增加2个中子；如果落在原子弹外部，就消失不见了，减少1个中子。
\\
\vspace{0.3cm}
一轮之后，我们统计中子总数，除以1000，如果这个比例大于1，就意味着下一轮核反应后中子数会继续增大，中子数呈幂级数增大，原子弹很快会爆炸。
\\
\vspace{0.3cm}
所以比例为1时，对应的就是球形原子弹的临界尺寸，单位是平均自由程。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>

#include"library.h"

const int n = 1000;

double ratio( double R, int n ){
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{链式反应：原子弹临界尺寸}
\lstset{firstnumber=last}
\begin{lstlisting}        
        int i, count = n;
        double r, costheta, phi, x, y, z;

        srand(time(NULL));
        for(i=0;i<n;i++){
                //1. 取初始位置
                x = R; y = R; z = R;//这些值没有什么意义，只是方便进入下面的while循环
                while( x*x + y*y + z*z > R*R ){
                        x = (double)rand()/RAND_MAX * 2 * R - R; // [-R,R]内均匀分布的随机数
                        y = (double)rand()/RAND_MAX * 2 * R - R;
                        z = (double)rand()/RAND_MAX * 2 * R - R;
                        //取得的点在边长为2R的立方体内均匀分布
                        //如果与原点距离大于R，会再次循环，重新取点
                        //如果与原点距离小于R，则结束while循环，得到球体内均匀分布的随机点
                }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{链式反应：原子弹临界尺寸}
\lstset{firstnumber=last}
\begin{lstlisting}                 
                //2. 自由飞行
                r = (double)rand()/RAND_MAX * 2;//飞行距离为[0,2]内均匀分布的随机数
                // costheta在[-1,1]内均匀分布，phi在[0,2pi]均匀分布，这样取点，可以得到均匀分布的立体角
                costheta = (double)rand()/RAND_MAX * 2 -1;
                phi = (double)rand()/RAND_MAX * 2 * M_PI;
                // (x,y,z) + \vec{r} = 终点
                x += r * sqrt( 1 - costheta * costheta ) * cos(phi);
                y += r * sqrt( 1 - costheta * costheta ) * sin(phi);
                z += r * costheta;
                //3. 判断终点位置
                if( x*x + y*y + z*z <= R*R ){
                // 终点在球内，触发一次反应，得到2个自由中子，即自由中子数+1
                        count ++;
                }
                else{
                // 终点在球外，消失不见，核弹内中子数-1
                        count --;
                }      
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{链式反应：原子弹临界尺寸}
\lstset{firstnumber=last}
\begin{lstlisting} 
        }
        return (double)count/n;
} 

double f(double R){// 这个函数的零点，对应 ratio=1
        return ratio(R, n) - 1;
}

int main(){
        
        double Rc = bisection(f, 0, 10, 1E-6);//二分法，在R=[0,10]区间内寻找临界半径
        cout<<"The critical radius R = "<<Rc<<" (unit is mean free path of neutron)"<<endl;
        return 0;
}     
\end{lstlisting}
可以得到，临界半径大约是1.36个平均自由程。
\end{frame}

\begin{frame}
	\begin{itemize}
		\item {\color{blue} 第十二章：蒙特卡洛方法 }
            \begin{itemize}
            \item {\color{blue} 简介}
            \item {\color{blue}蒙特卡洛模拟示例}
            \item {\color{blue}蒙特卡洛积分}
            \item Markov链
            \end{itemize}
		\item 第十三章：变分蒙特卡洛
	\end{itemize}
\end{frame}

\begin{frame}[fragile]{置信度与置信区间}
如果一个随机数服从正态分布$N(\mu, \sigma^2)$，那么它的概率密度函数是
\begin{equation}
f(x) = \frac{1}{\sqrt{2\pi}\sigma} exp\left(-\frac{(x-\mu)^2}{2\sigma^2} \right).
\end{equation}
那么，随机抽样落在$[\mu-\sigma, \mu+\sigma]$内的概率为
\begin{eqnarray}
F(\sigma) &=& \int^{\mu+\sigma}_{\mu-\sigma} \frac{1}{\sqrt{2\pi}\sigma} exp\left(-\frac{(x-\mu)^2}{2\sigma^2} \right) dx \\
&=& \int^{1}_{-1} \frac{1}{\sqrt{2\pi}} exp\left(-\frac{(x)^2}{2} \right) dx \\
&=& 68.268949 \%
\end{eqnarray}
类似地，随机抽样落在平均数左右2个标准差内的概率为$95.449974\%$，落在左右3个标准差内的概率为$99.730020\%$，落在左右4个标准差内的概率为$99.993666\%$。
\end{frame}

\begin{frame}[fragile]{置信度与置信区间}
在蒙特卡洛方法中，我们的计算结果$X$可能是随机的，服从分布
\begin{equation}
X \sim N(\mu, \sigma^2)
\end{equation}
可能理论上我们知道$\mu$就是物理量的理论值，但我们不知道它是多少，我们手上只有$n$次抽样的结果$x_i$，根据中心极限定理，我们知道抽样平均$\bar{x}$服从分布$N(\mu, \sigma^2/N)$，那么我们知道，抽样平均落在$[\mu - \sigma/\sqrt{n}, \mu + \sigma/\sqrt{n}]$的概率是68.268949\%，而
\begin{equation}
\mu - \sigma/\sqrt{n} \leq \bar{x} \leq \mu + \sigma/\sqrt{n}
\end{equation}
等价于
\begin{equation}
\bar{x} - \sigma/\sqrt{n} \leq \mu \leq \bar{x} + \sigma/\sqrt{n},
\end{equation}
所以我们给出的结果是：
\begin{shaded*}
{\color{blue}$\mu \in [\bar{x} - \sigma/\sqrt{n}, \bar{x} + \sigma/\sqrt{n}]$的概率是68.268949\%}。
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{置信度与置信区间}
68.268949\%的可能性，似乎不太靠谱，所以我们也可以改说：
\begin{shaded*}
{\color{blue}$\mu \in [\bar{x} - 4\sigma/\sqrt{n}, \bar{x} + 4\sigma/\sqrt{n}]$的概率是99.993666\%}。
\end{shaded*}
99.993666\%看起来挺可靠了，所以我们再努努力，增大抽样次数$n$，把区间缩小点，这个结论就很强大了。
\\
\vspace{0.3cm}
$[\bar{x} - 4\sigma/\sqrt{n}, \bar{x} + 4\sigma/\sqrt{n}]$：{\bf 置信区间}
\\
\vspace{0.3cm}
相应的概率99.993666\%：{\bf 置信度}。
\end{frame}

\begin{frame}[fragile]{简单蒙卡积分}
如果我们要求数值积分
\begin{equation}
I_f = \int^b_a f(x) dx,
\end{equation}
产生$[a,b]$内均匀分布的随机数，做$N$次抽样得到$x_i$，然后计算$f_i \equiv f(x_i)$，再求平均。
可以想象，如果随机数足够多，那么$[a,b]$内几乎所有小区域都被随机数均匀地覆盖到，那么应有
\begin{equation}
\frac{1}{N} \sum_i f_i \approx \frac{1}{b-a}I_f,
\end{equation}
所以有
\begin{equation}
I_f \approx (b-a) \bar{f},
\end{equation}
其中$\bar{f} = \frac{1}{N} \sum_i f_i$。
\end{frame}

\begin{frame}[fragile]{简单蒙卡积分}
你也可以把$f_i$看做是一个随机抽样，那么根据中心极限定理，$N$比较大的时候，$\bar{f}$服从$N(E(f), D(f)/N)$正态分布。

$f_i$的期望值是
\begin{equation}
E(f) = \frac{1}{b-a} \int^b_a f(x) dx,
\end{equation}
那么它的方差是
\begin{equation}
D(f) = \frac{1}{b-a} \int^b_a (f(x)-E(f))^2 dx = \frac{1}{b-a} \int^b_a f(x)^2 dx - E(x)^2.
\end{equation}
我们常常用抽样结果得到$D(f)$的近似估计，即
\begin{equation}
D(f) = \frac{1}{b-a} \int^b_a f(x)^2 dx - E(x)^2 \approx \overline{f^2} - \overline{f}^2.
\end{equation}
\end{frame}

\begin{frame}[fragile]{简单蒙卡积分}
这样，我们就可以用前面介绍的置信区间、置信度那一套办法，得到对积分的估计：
\begin{shaded*}
{\color{blue} $I_f \in (b-a)[\bar{f} - 4 \sqrt{ \frac{\overline{f^2} - \overline{f}^2}{n} }, \bar{f} + 4 \sqrt{ \frac{\overline{f^2} - \overline{f}^2}{n} }]$的可能性是99.993666\%。}
\end{shaded*}
这样，我们就得到了$I_f$的估计。
不出太大的意外的话，它都落在上面的区间内，我们可以多进行一些抽样，把区间限定得小一点，这个估计就很有用了。
\end{frame}

\begin{frame}[fragile]{例：求单位半圆面积}
单位半圆面积为
\begin{equation}
\int^1_{-1} \sqrt{1-x^2} dx = \pi/2,
\end{equation}
如果用蒙特卡洛积分来做，则有如下代码
\begin{lstlisting}
/*
 * 蒙特卡洛积分：求单位半圆面积
 */

#include<iostream>
using namespace std;
#include<cmath>

int main(){

        int i, n=1E8;
        double x, f, sum_f = 0, sum_f2 = 0, ave_f, ave_f2, I1, I2;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{例：求单位半圆面积}
\lstset{firstnumber = last}
\begin{lstlisting}
        for(int i=0;i<n;i++){
                x = (double)rand()/RAND_MAX*2 - 1;//取随机数 x 在 [-1,1] 均匀分布
                f = sqrt(1-x*x);// f*f + x*x =1
                sum_f += f;
                sum_f2 += f*f;
        }
        ave_f = sum_f / n;
        ave_f2 = sum_f2 / n;

        I1 = 2 * ( ave_f - 4 * sqrt( (ave_f2 - ave_f * ave_f)/n ) );//置信区间左边界
        I2 = 2 * ( ave_f + 4 * sqrt( (ave_f2 - ave_f * ave_f)/n ) );//置信区间右边界
        cout<<" If is in ["<<I1<<","<<I2<<"] by "<<99.993666<<"% chance."<<endl;
        // 4 sigma 置信度为 99.993666%
        cout<<" Exact value of If = "<<M_PI/2<<endl;
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{多重蒙特卡洛积分}
如果是多重积分
\begin{equation}
I_f = \int_V f dV,
\end{equation}
我们也可以用蒙特卡洛积分的方法，上文中的公式几乎是不变的，
\begin{shaded*}
{\color{blue} $I_f \in V[\overline{f} - 4 \sqrt{ \frac{\overline{f^2} - \overline{f}^2}{n} }, \bar{f} + 4 \sqrt{ \frac{\overline{f^2} - \overline{f}^2}{n} }]$的可能性是99.993666\%。}
\end{shaded*}
这个区间的长度正比于$\sqrt{1/n}$，与积分重数没有关系，这给蒙特卡洛方法带来优势。
这意味着，只要抽样次数一致，积分重数增长时，积分精度不变。
而在传统的积分方法中，要得到同样的精度，取点数按一维情况下的$m$次幂增长，$m$为积分重数。
\end{frame}

\begin{frame}[fragile]{重要抽样}
对于振荡剧烈的被积函数，如果在积分区域内均匀抽样，$\sqrt{ \overline{f^2} - \overline{f}^2 }$ 的值可能很大，导致置信区间很大，也就是精度不高。
我们可以把原来的积分稍作变形
\begin{equation}
I_f = \int_V f dV = \int_V \frac{f}{g} gdV,
\end{equation}
其中$g$是位置的函数。
这可以这么理解：
\begin{itemize}
\item[1]将每个体积微元扩大为原来的$g$倍，将该处的函数值$f$缩小为$f/g$，那么整体的积分不变。
\item[2]然后，我们可以进行随机抽样，得到$(f/g)_i$的值。
\item[3]我们在$gdV$这样的体积微元中抽样，相当于在$dV$中以概率密度函数$g$进行抽样，即在$V$中有的地方取点多一点，有的地方取点少一点。
\end{itemize}
\end{frame}

\begin{frame}[fragile]{重要抽样}
那么，相应地，我们得到积分的数值估计为
\begin{shaded*}
{\color{blue} $I_f \in \int_V g dV[\overline{f/g} - 4 \sqrt{ \frac{\overline{(f/g)^2} - \overline{(f/g)}^2}{n} }, \overline{f/g} + 4 \sqrt{ \frac{\overline{(f/g)^2} - \overline{(f/g)}^2}{n} }]$的可能性是99.993666\%。}
\end{shaded*}
如果$g$是概率密度函数，就有$\int_V g dV = 1$，当然，一般$V$比较规则，$g$比较解析，所以$\int_V g dV = 1$是容易保证的。则有
\begin{shaded*}
{\color{blue} $I_f \in [\overline{f/g} - 4 \sqrt{ \frac{\overline{(f/g)^2} - \overline{(f/g)}^2}{n} }, \overline{f/g} + 4 \sqrt{ \frac{\overline{(f/g)^2} - \overline{(f/g)}^2}{n} }]$的可能性是99.993666\%。}
\end{shaded*}
如果我们选取$g$，使得$f/g$的方差比较小（实际上，就是$g$尽量接近$f$的形状)，就会得到精度比较高的结果。
\end{frame}

\begin{frame}[fragile]{例：求单位半圆面积}
现在，我们实践重要抽样。选用
\begin{equation}
g(x) = \frac{2}{\pi} \sqrt{1-x*x},
\end{equation}
显然$g(x)$满足归一化条件
\begin{equation}
\int^1_{-1} g(x) dx = 1,
\end{equation}
而
\begin{equation}
f(x) = \sqrt{1-x*x},
\end{equation}
所以有
\begin{equation}
f(x) / g(x) = \frac{\pi}{2}
\end{equation}
\end{frame}

\begin{frame}[fragile]{例：求单位半圆面积}
那么，我们不用计算，就知道：
\begin{eqnarray}
\overline{f} &=& \pi/2, \\
\overline{f^2} &=& \pi^2/4 = \overline{f}^2.
\end{eqnarray}
带入前面的公式，我们知道：
\begin{shaded*}
{\color{blue} $I_f \in [\pi/2 - 0, \pi/2 + 0]$的可能性是99.993666\%。}
\end{shaded*}
实际上，如果你不用$4\sigma$的置信度，用更高的置信度，也会得到同样的置信区间，所以。。。
\begin{shaded*}
{\color{blue} $I_f \in [\pi/2 - 0, \pi/2 + 0]$的可能性是1。}
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{例：求单位半圆面积}
无论取样点数$n$等于多少，哪怕只是1，也会得到
\begin{shaded*}
{\color{blue} $I_f \in [\pi/2 - 0, \pi/2 + 0]$的可能性是1。}
\end{shaded*}

\begin{figure}
\includegraphics[width=0.6\textwidth]{why}
\end{figure}
\end{frame}

\begin{frame}[fragile]{例：求单位半圆面积}
这是因为，$g(x)$与$f(x)$完全相等，而$g(x)$已经归一化，意味着我们已经知道了积分的值。\\
\vspace{0.5cm}
实际应用的时候，我们是不知道积分的值的，所以无法使用这样的$g(x)$。\\
\vspace{0.5cm}
实际应用的时候，我们往往是用一个{\bf 形状非常接近 $f(x)$}、解析形式非常清楚的$g(x)$。\\
\vspace{0.5cm}
$f(x)$的形状越接近$g(x)$，抽样效率就越高。
\end{frame}

\begin{frame}[fragile]{例：重要抽样}
\begin{equation}
I_f = \int^{50}_0 (3-0.1*x)e^{-x} dx = 2.9,
\end{equation}
这个积分可以用普通抽样，也可以用重要抽样，下面我们展示同样的抽样次数对应的不同精度
\begin{lstlisting}
/*
 * 蒙特卡洛重要抽样积分： \int^50_0 (3-x)*exp(-x) dx
 */

#include<iostream>
using namespace std;
#include<cmath>

#include"library.h"

double f(double x){
        return (3-0.1*x)*exp(-x);
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{重要抽样}
\lstset{firstnumber = last}
\begin{lstlisting}
double g(double x){
        return exp(-x);
}

/*
 * 蒙特卡洛积分：简单抽样
 */
void simple_sampling(double xa, double xb, double (*f)(double x), int n){
        
        int i; 
        double x, y, sum_y = 0, sum_y2 = 0, ave_y, ave_y2, I1, I2;        
        for(int i=0;i<n;i++){
                x = (double)rand()/RAND_MAX*(xb-xa) + xa;//取随机数 x 在 [xa,xb] 均匀分布
                y = f(x);
                sum_y += y;
                sum_y2 += y*y;
        }
        ave_y = sum_y / n;
        ave_y2 = sum_y2 / n;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{重要抽样}
\lstset{firstnumber = last}
\begin{lstlisting}        
        I1 = (xb - xa) * ( ave_y - 4 * sqrt( (ave_y2 - ave_y * ave_y)/n ) );//置信区间左边界
        I2 = (xb - xa) * ( ave_y + 4 * sqrt( (ave_y2 - ave_y * ave_y)/n ) );//置信区间右边界
        cout<<"Simple sampling: integration value is "<< (I2+I1)/2<<" +- "<<(I2-I1)/2<<" by "<<99.993666<<"% chance."<<endl;
        // 4 sigma 置信度为 99.993666%
}

/*
 * 蒙特卡洛积分：重要抽样，抽样概率密度为 g(x)
 */
void importance_sampling(double xa, double xb, double (*f)(double x), double (*g)(double x), double gmax, int n){
        
        int i; 
        double x, y, sum_y = 0, sum_y2 = 0, ave_y, ave_y2, I1, I2;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{重要抽样}
\lstset{firstnumber = last}
\begin{lstlisting}        
        for(int i=0;i<n;i++){
                x = gnr_rand(xa, xb, g, gmax);
                y = f(x)/g(x);
                sum_y += y;
                sum_y2 += y*y;
        }
        ave_y = sum_y / n;
        ave_y2 = sum_y2 / n;

        I1 =   ave_y - 4 * sqrt( (ave_y2 - ave_y * ave_y)/n ) ;//置信区间左边界
        I2 =  ave_y + 4 * sqrt( (ave_y2 - ave_y * ave_y)/n ) ;//置信区间右边界
        cout<<"Importance sampling: integration value is "<< (I2+I1)/2<<" +- "<<(I2-I1)/2<<" by "<<99.993666<<"% chance."<<endl;
        // 4 sigma 置信度为 99.993666%
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{重要抽样}
\lstset{firstnumber = last}
\begin{lstlisting}
int main(){

        int n=1E6;
        simple_sampling(0, 50, f, n);
        importance_sampling(0, 50, f, g, 1, n);
        cout<<"exact value = 2"<<endl;

        return 0;
}
        
\end{lstlisting}
取$n=1E6$时得到结果如下
\begin{lstlisting}
Simple sampling: integration value is 2.89711 +- 0.0577714 by 99.9937% chance.
Importance sampling: integration value is 2.90012 +- 0.000399186 by 99.9937% chance.
exact value = 2.9
\end{lstlisting}
可见，重要抽样得到的结果更精确。
\end{frame}

\begin{frame}
	\begin{itemize}
		\item {\color{blue} 第十二章：蒙特卡洛方法 }
            \begin{itemize}
            \item {\color{blue} 简介}
            \item {\color{blue} 蒙特卡洛模拟示例 }
            \item {\color{blue} 蒙特卡洛积分 }
            \item {\color{blue} Markov链 }
            \end{itemize}
		\item 第十三章：变分蒙特卡洛
	\end{itemize}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链}
根据统计系综理论，正则系综（恒温系统）力学量的宏观测量值，等于不同微观状态做如下统计平均得到的结果，
\begin{eqnarray}
\langle A \rangle = \frac{ \int A e^{-\beta E} d \Omega }{ \int e^{-\beta E} d \Omega},
\end{eqnarray}
其中 $E$ 为微观状态对应的体系总能量，$\beta = 1/(k_B T)$，$k_B$为玻尔兹曼常数。\\
\vspace{0.3cm}
如果我们在整个正则系综里，以$e^{-\beta E}$的相对概率进行抽样，则有
\begin{eqnarray}
\langle A \rangle = \frac{1}{n} \sum_i A_i,
\end{eqnarray}
我们的计算机这么小，当然不可能抽取所有的微观状态。
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链}
我们能做的是构造一个有限的空间，其中有一些微观状态，它们的{\bf 相对概率满足正则分布}。
那么则有
\begin{eqnarray}
\langle A \rangle \approx \frac{1}{n} \sum_i A_i.
\end{eqnarray}

如果这些微观状态呈一个序列，这个序列就叫做{\bf 马尔科夫（Markov）链}。
构造这样的一个链以后，从链上抽样做平均，即可近似地计算体系的宏观物理量。
这样的模拟就叫做{\bf 麦氏模拟(Metropolis)}。
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链}
我们可以先做一个最简单的例子，来熟悉上面说的这些概念。\\
\vspace{0.3cm}
假设只有一个原子，自旋可以向上，也可以向下，外加一个向上的磁场，则体系哈密顿量为
\begin{equation}
H = - \mu B S_z,
\end{equation}
其中$\mu$为原子磁矩，所以自旋向上时的自旋、能量、相对概率为
\begin{equation}
S_z = 1, ~~ E_{up} = -\mu B, ~~e^{-E_{up}/(kT)} = e^{\mu B/(kT)},
\end{equation}
自旋向下的自旋、能量、相对概率为
\begin{equation}
S_z = -1, ~~ E_{down} = \mu B, ~~e^{-E_{down}/(kT)} = e^{-\mu B/(kT)}.
\end{equation}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链}
体系能量的测量值为
\begin{eqnarray}
\langle E \rangle = \frac{ \int E e^{-\beta E} d \Omega }{ \int e^{-\beta E} d \Omega}
= \frac{-\mu B e^{\mu B/(kT)} + \mu B e^{-\mu B/(kT)}}{e^{\mu B/(kT)} + e^{-\mu B/(kT)}},
\end{eqnarray}
体系自旋的测量值为
\begin{eqnarray}
\langle S_z \rangle = \frac{ \int S_z e^{-\beta E} d \Omega }{ \int e^{-\beta E} d \Omega}
= \frac{ e^{\mu B/(kT)} -  e^{-\mu B/(kT)}}{e^{\mu B/(kT)} + e^{-\mu B/(kT)}},
\end{eqnarray}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链}
不妨设$\mu B = 2$，则$\langle E \rangle$与$kT$的关系为

\begin{figure}
\includegraphics[width=0.6\textwidth]{src/mc_markov_demon_E}
\end{figure}

温度极低时，自旋向上，能量为$-2$；温度极高时，自旋一半概率向上，一半概率向下，能量平均为$0$。
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链}
$\langle S_z \rangle$与$kT$的关系为

\begin{figure}
\includegraphics[width=0.6\textwidth]{src/mc_markov_demon_Sz}
\end{figure}
温度极低时，自旋向上，$S_z=1$；温度极高时，自旋一半概率向上，一半概率向下，$S_z$自旋平均为$0$。
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链}
那么，给定温度，我们要构造的马尔科夫链如下
\begin{equation}
\uparrow \downarrow \uparrow \uparrow \downarrow \uparrow \uparrow \downarrow \uparrow \cdots
\end{equation}
上面的$\uparrow$出现的相对概率为$e^{\mu B/(kT)}$，$\downarrow$出现的相对概率为$-e^{\mu B/(kT)}$。
\\
\vspace{0.3cm}
然后，我们在链上随机抽样，做平均，就可以得到
\begin{eqnarray}
\langle E \rangle &\approx& \overline{E}, \\
\langle S_z \rangle &\approx& \overline{S_z}.
\end{eqnarray}
所以关键是用什么算法实现马尔科夫链
\end{frame}

\begin{frame}[fragile]{Metropolis算法}
如果链上某个点是微观状态$X$，则计算$E_{X'}, E_X$，按如下概率转换，
\begin{eqnarray}
T_{X \rightarrow X'} = \left\{
\begin{aligned}
&e^{ - \beta (E_{X'} - E_{X}) }, &E_{X'} > E_{X};\\
&1, ~~~~~~~~~~~~~~~~~~~~~~~~&E_{X'} < E_{X}.
\end{aligned}
\right.
\end{eqnarray}
可以看到
\begin{equation}
\frac{ T_{X' \rightarrow X} }{T_{X \rightarrow X'}} = e^{-\beta(E_X - E_{X'})},
\end{equation}
最终链上各种微观状态达到平衡以后，会有
\begin{equation}
\rho(X',t)T_{X' \rightarrow X} = \rho(X,t)T_{X \rightarrow X'},
\end{equation}
其中$\rho(X,t)$是在第$t$个节点微观状态$X$出现的概率。所以有
\begin{equation}
\frac{\rho(X,t)}{\rho(X',t)} 
= \frac{ T_{X' \rightarrow X} }{T_{X \rightarrow X‘}}
= \frac{e^{-\beta E_X}}{e^{-\beta E_{X'}}},
\end{equation}
\end{frame}

\begin{frame}[fragile]{Metropolis算法}
假设初始状态为$\uparrow$，因为$E_{down} > E_{up}$，所以
\begin{equation}
T_{\uparrow \rightarrow \downarrow} = e^{-\beta(E_{down} - E_{up})} = e^{-2\beta\mu B},
\end{equation}
于是我们取$[0,1]$均匀分布的随机数$\xi$；如果$\xi<e^{-2\beta\mu B}$，就翻转原子，得到下一个节点$\downarrow$；如果$\xi \geq e^{-2\beta\mu B}$，就不翻转，下一个节点仍为$\uparrow$。\\
\vspace{0.5cm}
如果节点为$\downarrow$，因为$E_{up} < E_{down}$，所以
\begin{equation}
T_{\downarrow \rightarrow \uparrow} = 1,
\end{equation}
于是我们翻转原子，得到下一个节点$\uparrow$。\\
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链}
\begin{lstlisting}
/*
 * 单原子自旋可以向上、向下，磁矩为mu，外场B向上，哈密顿量为
 *  H = - mu B Sz
 *  下面的代码输出马尔科夫链在屏幕上
 */

#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>

void markov_single_atom(int n, double muB, double kT, double & ave_E, double & ave_Sz){
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链}
\lstset{firstnumber = last}
\begin{lstlisting}
        double sum_E = 0, sum_Sz = 0;
        int Sz = 1;
        for(int i=0;i<n;i++){
                if(Sz ==1){//如果自旋向上
                        //cout<<"+";
                        sum_E += - muB;
                        sum_Sz += Sz;
                        if( (double)rand()/RAND_MAX  < exp(-2*muB/kT) )
                                Sz = -1;
                        //以exp(-2*muB/kT)的概率翻转为向下
                        //以1-exp(-2*muB/kT)的概率保持向上                      
                }
                else{//如果自旋向下
                        sum_E += muB;
                        sum_Sz += Sz;
                        //cout<<"-";
                        Sz = 1;
                        //翻转为向上
                }
        }
        cout<<endl;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链}
\lstset{firstnumber = last}
\begin{lstlisting}
        ave_E = sum_E /n;
        ave_Sz = sum_Sz /n;
}

int main(){

        int n=10000;
        double muB = 2;
        double kT, ave_E, ave_Sz;
        ofstream fp("mc_markov_demon_sampling.txt");
        for(kT = 0; kT <50; kT += 0.1){
                markov_single_atom(n, muB, kT, ave_E, ave_Sz);
                fp<< kT<<"\t"<<ave_E<<"\t"<<ave_Sz<<endl;
        }
        fp.close();
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链}
仍然设$\mu B = 2$，上面的代码运行出来得到$\langle E \rangle$与$kT$的关系为

\begin{figure}
\includegraphics[width=0.6\textwidth]{src/mc_markov_demon_sampling_E}
\end{figure}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链}
$\langle S_z \rangle$与$kT$的关系为
\begin{figure}
\includegraphics[width=0.6\textwidth]{src/mc_markov_demon_sampling_Sz}
\end{figure}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链：多个孤立原子}
下面我们再往前走一步，选择256个相同的孤立原子，都处在外场B中，
\begin{equation}
H = - \mu B \sum^{256}_{i=1}  s^{(i)}_z = - \mu B S_z,
\end{equation}
其中$s^{(i)}_z$表示第i个原子的自旋第3分量，$S_z$表示总自旋的第3分量。\\
\vspace{0.2cm}
用如下算法实现metropolis
\begin{shaded*}
\begin{itemize}
\item[1] 依次征询256个原子，依如下概率决定是否翻转。
\begin{eqnarray}
T_{X \rightarrow X'} = \left\{
\begin{aligned}
&e^{ - \beta (E_{X'} - E_{X}) }, &E_{X'} > E_{X};\\
&1, ~~~~~~~~~~~~~~~~~~~~~~~~&E_{X'} < E_{X}.
\end{aligned}
\right.
\end{eqnarray}
\item[2] 每征询完一轮256个原子，进行一次抽样。
\end{itemize}
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链：多个孤立原子}
只需对原来的程序略作修改
\begin{lstlisting}
/*
 * 单原子自旋可以向上、向下，磁矩为mu，外场B向上，哈密顿量为
 *  H = - mu B Sz
 *  下面的代码输出马尔科夫链在屏幕上
 */
#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>

void markov_256atom(int n, double muB, double kT, double & ave_e, double & ave_sz){

        double sum_E = 0, sum_Sz = 0;//大写的E、Sz表示体系能量、总自旋第3分量
        int sz[256]; for(int i=0;i<256;i++) sz[i] = 1;
        // 小写的 ave_e, ave_sz 表示平均每原子能量、自旋第3分量的系综平均
        double E, Sz;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链：多个孤立原子}
\lstset{firstnumber = last}
\begin{lstlisting}
        for(int i=0;i<n;i++){
                Sz=0; for(int j=0;j<256;j++) Sz += sz[j];
                sum_Sz += Sz;
                sum_E += - muB * Sz;
                for(int j=0;j<256;j++){
                        if(sz[j] ==1){
                                if( (double)rand()/RAND_MAX  < exp(-2*muB/kT) )
                                        sz[j] = -1;
                        }
                        else{
                                sz[j] = 1;
                        }
                }
        }

        ave_e = sum_E /n/256;
        ave_sz = sum_Sz /n/256;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链：多个孤立原子}
\lstset{firstnumber = last}
\begin{lstlisting}
int main(){

        int n=1000;
        double muB = 2;
        double kT, ave_e, ave_sz;
        ofstream fp("mc_markov_demon_sampling_256atom.txt");
        for(kT = 0; kT <50; kT += 0.1){
                markov_256atom(n, muB, kT, ave_e, ave_sz);
                fp<< kT<<"\t"<<ave_e<<"\t"<<ave_sz<<endl;
        }
        fp.close();
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链：多个孤立原子}
结果当然和单原子的结果是一样的
\begin{figure}
\includegraphics[width=0.6\textwidth]{src/mc_markov_demon_256atom_E}
\end{figure}
\end{frame}

\begin{frame}[fragile]{马尔科夫(Markov)链：多个孤立原子}
$\langle s_z \rangle$与$kT$的关系为
\begin{figure}
\includegraphics[width=0.6\textwidth]{src/mc_markov_demon_256atom_Sz}
\end{figure}
\end{frame}

\begin{frame}[fragile]{2维Ising模型}
现在，我们将256个原子按照 16 x 16 排成方阵
\begin{figure}
\includegraphics[width=0.7\textwidth]{src/IsingModel_demon}
\end{figure}
\end{frame}

\begin{frame}[fragile]{2维Ising模型}
不仅考虑外加磁场，并且考虑相邻原子之间的相互作用
\begin{equation}
H = - J \sum_{\{i,j\}} s^{(i)}_z s^{(j)}_z - \mu B \sum_i s^{(i)}_z,
\end{equation}
其中，$\{i,j\}$表示$i,j$相邻。\\
\vspace{0.3cm}
如果$J>0$，那么相邻原子取向一致时，体系能量更低；如果$J<0$，那么相邻原子取向相反时，体系能量更低。\\
\vspace{0.3cm}
\end{frame}

\begin{frame}[fragile]{2维Ising模型：周期性边界条件}
我们定义了边长$L=16$，将原子标上序号$0,1,\cdots,L*L-1$，则第$i$个原子的上邻居为
\begin{lstlisting}
i>(L-1)*L-1 ? i-(L-1)*L : i+L
\end{lstlisting}
下邻居为
\begin{lstlisting}
i<L? i+(L-1)*L : i-L
\end{lstlisting}
左邻居为
\begin{lstlisting}
i%L ==0 ? i+L-1 : i-1
\end{lstlisting}
右邻居为
\begin{lstlisting}
(i+1)%L==0 ? i-(L-1) : i+1
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维Ising模型：翻转原子}
那么，第$i$个原子的所有邻居的$s_z$总和为
\begin{lstlisting}
S_neighbour = spin[ (i>(L-1)*L-1 ? i-(L-1)*L : i+L) ] + spin[ (i<L? i+(L-1)*L : i-L) ]
              + spin[ (i%L ==0 ? i+L-1 : i-1) ] + spin[ ((i+1)%L==0 ? i-(L-1) : i+1) ];
\end{lstlisting}
翻转第$i$个原子以后，体系能量变化为
\begin{lstlisting}
dE= 2 * J * S_neighbour * spin[i] + 2 * muB * spin[i];
\end{lstlisting}
\begin{shaded*}
\begin{itemize}
\item[1] 依次征询256个原子，依如下概率决定是否翻转。
\begin{eqnarray}
T_{X \rightarrow X'} = \left\{
\begin{aligned}
&e^{ - \beta (E_{X'} - E_{X}) }, &E_{X'} > E_{X};\\
&1, ~~~~~~~~~~~~~~~~~~~~~~~~&E_{X'} < E_{X}.
\end{aligned}
\right.
\end{eqnarray}
\item[2] 每征询完一轮256个原子，进行一次抽样。
\end{itemize}
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{2维Ising模型}
\begin{lstlisting}
/*
 * 2-dimensional Ising model, with Markov chain in Monte Carlo method
 *      H = - J \sum_{ij} sz_i sz_j - muB \sum_i sz_i,
 * where J is the interaction strength between two atoms, negative/positive for ferromagnetic/anti-ferromagnetic material, k is the Boltzman constant, T is the temperature, we set J as unit.
 */

#include<iostream>
#include<stdlib.h>
#include<time.h>
#include<fstream>
#include<cmath>
using namespace std;

void Ising2d(int length, int L, double kT, double muB, double & sz_ave, double & e_ave, double & e2_ave){

        long i,j,k, Sz, S_neighbour, count=0;
        double J=1, E, dE, sz_sum=0, e_sum=0, e2_sum=0;
        int *spin=new int [L*L];//spin of each atom, which can only be 1 or -1
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维Ising模型}
\lstset{firstnumber = last}
\begin{lstlisting}
        bool flag_flip;
        for(i=0;i<L*L;i++) //initially the lattice is magnetizied.
                spin[i]=1;

        srand((unsigned)time(0));//initialize rand() with a random seed, according to the real time when this code runs.
        Sz= L*L;//initial value of total spin
        E=-J*2*L*L - muB * Sz;//initial value of total energy of the system

        for(j=0;j<length;j++){

                for(i=0;i<L*L;i++){

                        S_neighbour = spin[ (i>(L-1)*L-1 ? i-(L-1)*L : i+L) ] + spin[ (i<L? i+(L-1)*L : i-L) ]
                                + spin[ (i%L ==0 ? i+L-1 : i-1) ] + spin[ ((i+1)%L==0 ? i-(L-1) : i+1) ];
                        dE= 2 * J * S_neighbour * spin[i] + 2 * muB * spin[i];
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维Ising模型}
\lstset{firstnumber = last}
\begin{lstlisting}
                        flag_flip=false;
                        if(dE<0){
                               flag_flip=true;
                        }
                        else{
                                if( (double)rand()/RAND_MAX < exp(-dE/kT) ){// rand()/RAND_MAX uniformly distributes in [0,1], so it has possibility exp(-dE/kT) to be less than exp(-dE/kT)
                                        flag_flip=true;
                                }
                        }

                        if(flag_flip){
                                E += dE;
                                spin[i] *= -1;//flip the spin
                                Sz += 2*spin[i];// total spin and total energy changes
                        }
                }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维Ising模型}
\lstset{firstnumber = last}
\begin{lstlisting}
                // after one round, collects data
                if(j>0.1*length){// We assume the system is in equilibrium after 0.1*length
                        sz_sum += (double)Sz/(L*L);
                        e_sum += E/(L*L);
                        e2_sum += E*E/(L*L);

                        count ++;
                }
        }

        sz_ave = sz_sum/count;
        e_ave = e_sum/count;
        e2_ave = e2_sum/count;

        delete [] spin;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维Ising模型}
\lstset{firstnumber = last}
\begin{lstlisting}
int main(){

        int length = 1E4, L=16;
        double sz_ave, e_ave, e2_ave, muB=0;

        ofstream fp("Ising2d.txt");
        for(double kT=0.1; kT<=5; kT+=0.1){
                Ising2d(length, L, kT, muB, sz_ave, e_ave, e2_ave);
                fp<< kT <<"\t"<< sz_ave << "\t" << e_ave <<endl;
        }
        fp.close();

        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{2维Ising模型}
另有一个相似的代码，输出一个gnuplot脚本，产生动图，进行演示。
下面是$\langle e \rangle$与$kT/J$之间的关系
\begin{figure}
\includegraphics[width=0.7\textwidth]{src/IsingModel_ave_e}
\end{figure}
\end{frame}

\begin{frame}[fragile]{2维Ising模型}
以及$\langle sz \rangle$与$kT/J$之间的关系
\begin{figure}
\includegraphics[width=0.7\textwidth]{src/IsingModel_ave_sz}
\end{figure}
可以看到，$kT/J$在$2.2$附近有一个剧烈的变化点，这个点对应着铁磁-非铁磁相变。
\end{frame}

\begin{frame}
	\begin{itemize}
		\item {\color{blue} 第十二章：蒙特卡洛方法 }
            \begin{itemize}
            \item {\color{blue} 简介}
            \item {\color{blue}蒙特卡洛模拟示例}
            \item {\color{blue}蒙特卡洛积分}
            \item {\color{blue}Markov链}
            \end{itemize}
		\item {\color{blue}第十三章：变分蒙特卡洛}
	\end{itemize}
\end{frame}

\begin{frame}[fragile]{变分法}
在初等量子力学中，我们学过变分法，用来优化体系的一个近似本征态，下面简单回顾这一理论。
若体系的哈密顿量 $\hat{H}$ 已知，存在本征值 $\epsilon_i$ 与 $\phi_i$，使得
\begin{eqnarray}
\hat{H} \phi_i = \epsilon_i \phi_i,~~ i = 1, 2, 3, \cdots.
\end{eqnarray}
如果问题复杂，无法严格解出 $\epsilon_i, \phi_i$，那么能得到基态能量和基态波函数的近似解，就已经很满足了。

于是我们构造一个波函数 $\Psi$。无论怎样构造，它一定可以表示成 $\phi_i$ 的线性展开（因为 $\hat{H}$ 的所有本征态构成完备基），
\begin{eqnarray}
\Psi = a_1 \phi_1 + a_2 \phi_2 + \cdots,
\end{eqnarray}
虽然我们不知道 $\phi_1, \phi_2, \cdots $ 长什么样子，$\Psi$ 是它们的线性展开，这一点总是没错的。
\end{frame}

\begin{frame}[fragile]{变分法}
另外，
\begin{eqnarray}
\langle E \rangle &=& \frac{ \langle \Psi | \hat{H} | \Psi \rangle }{ \langle \Psi | \Psi \rangle }
= \frac{ \sum\limits^n_{i,j=1} a^*_i a_j \langle \phi_i | \hat{H} | \phi_j \rangle }{ \sum_i |a_i|^2 }
= \frac{ \sum\limits^n_{i,j=1} a^*_i a_j \epsilon_j \delta_{ij} }{ \sum_i |a_i|^2 }
\nonumber \\
&=& \frac{ \sum\limits^n_{i=1} |a_i|^2 \epsilon_i }{ \sum_i |a_i|^2 }
\geq \frac{ \sum\limits^n_{i=1} |a_i|^2 \epsilon_1 }{ \sum_i |a_i|^2 }
= \epsilon_1
\end{eqnarray}
当且仅当 $a_2 = a_3 = \cdots =0$ 时，上式中的 $\geq$ 符号取等号。
也就是说，
\begin{itemize}
\item[1] 对于任意不为0的试探波函数，$ \langle E \rangle$ 一定不小于体系基态能量。
\item[2] $| \Psi \rangle$ 越接近体系基态波函数，$ \langle E \rangle$ 就会越接近体系基态能量。
\item[3] $| \Psi \rangle$ 往往有一个带参数的拟设形式，在此形式下，$ \langle E \rangle$ 往往是可以计算的，所以我们要做的，就是调节 $\Psi$（即调节其中参数），使得 $ \langle E \rangle$ 尽量小，取极小值的时候，就把 $\Psi$ 看作 $\phi_1$ 的近似，把 $ \langle E \rangle$ 看作是 $\epsilon_1$ 的近似。
\end{itemize}
\end{frame}

\begin{frame}[fragile]{变分蒙特卡洛(VMC)}
在给定参数$\alpha$下，我们可以直接做积分，然后计算
\begin{equation}
\langle E \rangle_\alpha = \frac{ \langle \Psi_\alpha | \hat{H} | \Psi_\alpha \rangle }{ \langle \Psi_\alpha | \Psi_\alpha \rangle },
\end{equation}
然后改变$\alpha$，使得上式达到极小值。
$\Psi_\alpha(\vec{r})$在$\vec{r}$空间上往往是变化很大的，在大部分区域接近于0，可以使用重要抽样。

我们也可以用Markov链的思路，换一种方法求取$\langle E \rangle_\alpha$。

首先，在$\vec{r}$区间内随机选择多个随机行走者（random walker）。
对于每一个 random walker，计算
\begin{equation}
E_L (\vec{r})= H \Psi (\vec{r}) / \Psi(\vec{r}),
\end{equation}
这个量叫做局域能量，我们可以写出
\begin{equation}
\langle E \rangle = \frac{ \int \Psi(\vec{r})^2 E_L(\vec{r}) d\vec{r} }{ \int \Psi(\vec{r})^2  d\vec{r} },
\end{equation}
所以，相当于以$\Psi(\vec{r})^2$为相对概率，对$E_L(\vec{r})$进行积分。
所以我们可以构造相对概率为$\Psi(\vec{r})^2$的马尔科夫链，对$E_L(\vec{r})$进行抽样，得到$\langle E \rangle$。
\end{frame}

\begin{frame}[fragile]{变分蒙特卡洛(VMC)}
\begin{shaded*}
\begin{itemize}
\item[1] 在$\vec{r}$空间均匀选取大量random walker。
\item[2] 每个 random walker 计算
\begin{equation}
\langle E \rangle_\alpha = \frac{ \langle \Psi_\alpha | \hat{H} | \Psi_\alpha \rangle }{ \langle \Psi_\alpha | \Psi_\alpha \rangle },
\end{equation}
\item[3] 每个 random walker 产生一个随机行走步长，比如$x,y,z$都有一般概率加上$\delta$，一般概率减去$\delta$；也可以做得更均匀一点，让步长为朝随机立体角行走$\delta$。
然后计算终点的波函数值　$\Psi(\vec{r}^\prime)$，并计算比例
\begin{equation}
p = \Psi(\vec{r}^\prime)^2 / \Psi(\vec{r})^2,
\end{equation}
\end{itemize}
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{变分蒙特卡洛(VMC)}
\begin{shaded*}
\begin{itemize}
\item[4] 根据如下概率，决定是否接受这次随机行走
\begin{eqnarray}
T_{\vec{r} \rightarrow \vec{r}^\prime} = \left\{
\begin{aligned}
&p, &p<1;\\
&1, &p>1.
\end{aligned}
\right.
\end{eqnarray}
\item[5]根据细致平衡条件，最终random walker出现在各处的概率与该处波函数平方成正比。
所以我们做$E_L$的抽样即可，
\begin{equation}
\langle E \rangle_\alpha \approx \overline{E_L}.
\end{equation}
\item[6]变化$\alpha$的值，重复上述过程，做出$\overline{E_L} \sim \alpha$的曲线，找到极小值点，即完成近似求解。
\end{itemize}
\end{shaded*}
\end{frame}

\begin{frame}[fragile]{VMC:一维谐振子}
没有量纲的一维谐振子薛定谔方程为
\begin{equation}
H \psi(x) = (-\frac{1}{2} \frac{d^2}{dx^2} + \frac{1}{2} x^2) \psi(x),
\end{equation}
基态的严格解为
\begin{equation}
\psi_{g.s.} (x) = A e^{-x^2/2},
\end{equation}
其中$A$为归一化因子。
基态能量为$e_{g.s.} = 1/2$。

为了演示VMC的思路，我们选择拟设形式
\begin{equation}
\psi_\alpha (x) = e^{-\alpha x^2},
\end{equation}
则局域能量为
\begin{equation}
E_L(x) = H \psi(x) / \psi(x) = \alpha + (\frac{1}{2} - 2\alpha^2) x^2,
\end{equation}
利用这个公式，如下代码实现VMC，
\end{frame}

\begin{frame}[fragile]{VMC:一维谐振子}
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>

double psi(double alpha, double x){
        return exp(-alpha*x*x);
}
double EL(double alpha, double x){
        return alpha + (0.5-2*alpha*alpha)*x*x;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{VMC:一维谐振子}
\lstset{firstnumber = last}
\begin{lstlisting}
/*
 * vmc1d_aveEL:         变分蒙特卡洛方法，返回<EL>值
 * doubla alpha:        参数
 * int n:               number of random walker
 * double (*psi)()      波函数 psi(alpha, x)
 * double (*EL)()       局域能量 H psi(alpha, x) / psi(alpha, x)
 * double a, b:         x　轴定义域, 所有random walker被限制在[a,b]内
 * double h:            随机行走步长
 * int N:               每个 random walker的markov链长度
 */
double vmc1d_aveEL(double alpha, int n,
                double (*psi)(double alpha, double x), double (*EL)(double alpha, double x),
                double a, double b, double h, int N){

        srand(time(NULL));
        int i,j,count;
        double x, d, y, z, p;
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{VMC:一维谐振子}
\lstset{firstnumber = last}
\begin{lstlisting}
        y=0;
        for(i=0;i<n;i++){
                x = (double)rand()/RAND_MAX * (b-a) + a;//初始位置
                z=0;
                count=0;
                for(j=0;j<N;j++){
                        if( j > 0.1 *N ){
                                z += EL(alpha, x);
                                count ++;
                        }
                        d = ((rand() % 2)*2-1)*h; // -h or h
                        if( x+d > b || x+d < a ) continue;
                        p = pow( psi(alpha, x+d), 2) / pow( psi(alpha,x), 2);
                        if(p>1) x += d;
                        else{
                                if( (double)rand()/RAND_MAX < p ) x+=d;
                        }
                }
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{VMC:一维谐振子}
\lstset{firstnumber = last}
\begin{lstlisting}
                }
                z /= count;
                y += z;
        }
        return y/n;
}

int main(){

        ofstream fp("vmc1dho_aveEL.txt");
        for(double alpha=0.4; alpha <=0.6; alpha += 0.01){
                fp << alpha << "\t" << vmc1dho_aveEL(alpha, 100, 0.2, 10000) <<endl;
        }
        fp.close();
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{VMC:一维谐振子}
vmc1d\_aveEL函数可以复制粘贴进课程的 library.h　文件，后面就可以直接调用这个函数，不用再次编写了。\\


代码运行以后，通过以下gnuplot代码可以生成图片
\begin{lstlisting}
set term png
set output "vmc1dho.png"
set xlabel "{/Symbol a}"
set ylabel "<E>"
plot "vmc1dho_aveEL.txt" u 1:2 w l lw 2 t "VMC: 1d Harmonic Oscillator"
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{VMC:一维谐振子}
得到下图。
可以看到，最低点对应着$\alpha=1/2$，即体系基态波函数，相应的$\langle E \rangle = 1/2$对应着基态能量$1/2$。
\begin{figure}
\centering
\includegraphics[width=0.6\textwidth]{src/vmc1dho}
\caption{变分蒙特卡洛：一维谐振子拟设波函数$\psi_\alpha(x) = e^{-\alpha x^2}$，在$\alpha$取$0.4\sim 0.6$时，得到的$<E> = <\psi_\alpha|H|\psi_\alpha>/<\psi_\alpha | \psi_\alpha>$值。}
\label{fig:vmc1dho}
\end{figure}
\end{frame}

\begin{frame}[fragile]{VMC:氢原子}
没有量纲的氢原子薛定谔方程为
\begin{equation}
(-\frac{1}{2}\nabla^2 - \frac{1}{r} ) \psi(\vec{r}) = E \psi(\vec{r}),
\end{equation}
分离变量以后，设
\begin{equation}
\psi(\vec{r}) = ( u(r)/r ) Y_{lm}(\theta, \phi),
\end{equation}
$s$波方程转化为一维径向方程
\begin{equation}
(-\frac{1}{2} \frac{d^2}{dr^2} - \frac{1}{r} ) u(r) = E u(r),
\end{equation}
基态的严格解为
\begin{equation}
u(r) = A r e^{-r},
\end{equation}
对应着
\begin{equation}
\psi_{g.s.} (\vec{r}) = A e^{-r}.
\end{equation}
其中$A$为归一化因子。
基态能量为$e_{g.s.} = -1/2$。
\end{frame}

\begin{frame}[fragile]{VMC:氢原子}
为了演示VMC的思路，我们选择拟设形式
\begin{equation}
u(r) = r e^{-\alpha r},
\end{equation}
没有归一化，但是不影响我们后面的计算，因为我们真正使用的是波函数在不同地方的相对大小。

局域能量为
\begin{equation}
E_L(r) = (-\frac{1}{2} \frac{d^2}{dr^2} - \frac{1}{r} )　u(r) / u(r) = -\frac{1}{2}\alpha^2 + (\alpha-1)/r,
\end{equation}
利用这个公式，如下代码实现VMC，
\end{frame}

\begin{frame}[fragile]{VMC:氢原子}
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>
#include"library.h"
double u(double alpha, double r){
        return r*exp(-alpha*r);
}
double EL(double alpha, double r){
        return -0.5*alpha*alpha + (alpha-1)/r;
}
int main(){

        ofstream fp("vmcHydrogen_aveEL.txt");
        for(double alpha=0.8; alpha <=1.2; alpha += 0.01){
                fp << alpha << "\t" << vmc1d_aveEL(alpha, 100, u, EL, 0, 5, 0.2, 10000) <<endl;
        }
        fp.close();
        return 0;
}
\end{lstlisting}
\end{frame}

\begin{frame}[fragile]{VMC:氢原子}
得到下图。
可以看到，最低点对应着$\alpha=1$，即体系基态波函数为$e^{-r}$，相应的$\langle E \rangle = -1/2$对应着基态能量$-1/2$。
\begin{figure}
\centering
\includegraphics[width=0.6\textwidth]{src/vmcHydrogen}
\caption{变分蒙特卡洛：氢原子拟设波函数$\psi_\alpha(\vec{r}) = e^{-\alpha r}$，在$\alpha$取$0.8\sim 1.2$时，得到的$<E> = <\psi_\alpha|H|\psi_\alpha>/<\psi_\alpha | \psi_\alpha>$值。}
\label{fig:vmcHydrogen}
\end{figure}
\end{frame}

\begin{frame}[fragile]{作业}
\begin{itemize}
\item [1] 用点名程序点一百万次名，并统计每个人被点到的次数，做图观察并给出自己的结论：是否每个人被点到的频率差不多。
\item [2] 产生一百万个$[0,1]$区间内的随机数，并统计它们分布于$[0,0.1],\cdots,[0.9,1]$区间内的频率，观察是否都差不多。
\item [3] 产生$[0,\infty]$内连续分布的随机数，概率密度函数为
\begin{equation}
f(x) = e^{-x}. \nonumber
\end{equation}
试在$[0,10]$内产生1百万个这样的随机数，并统计各处频率，与概率密度函数比照，检验代码正确性。
\item[4] 用蒙特卡洛积分计算$\int^{\pi}_0 \sin(x) dx$，给出$4\sigma$精度的置信区间。抽样次数自己决定。
\item[5] 对单原子Markov链的代码稍作修改，统计$\mu B =2, kT = 4$时，Markov链上自旋向上、向下的相对概率，并与理论结果$\rho(down)/\rho(up) = exp(-2\mu B / kT)$相比较。
\end{itemize}
\end{frame}

\begin{frame}[fragile]{作业}
\begin{itemize}
\item [*6] 自己编写3维Ising模型，并寻找相变点。
\item [7] 尝试用VMC方法，用试探波函数$\psi(x) = x e^{-\alpha x^2}$，求解1维势井基态：
\begin{equation}
V(x) = \left\{
\begin{aligned}
&x^2/2, ~~~ 0<x<\infty\\
&\infty, ~~ x<0
\end{aligned}
\right.
\end{equation}
\end{itemize}
薛定谔方程变为
\begin{equation}
(- \frac{d^2}{dx^2} + V(x) ) \psi = E \psi.
\end{equation}
理论解为
\begin{eqnarray}
E_{g.s.} &=& 1.5, \\
\psi_{g.s.} &\propto& x e^{-x^2/2}.
\end{eqnarray}
\end{frame}

\end{document}